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Abstract. The general PDE governing linear, adiabatic, nonraradial oscillations in a spherical, differentially and 
slowly rotating non-magnetic star is derived. This equation describes mainly low-frequency and high-degree g- 
modes, convective (j-modes, and rotational Rossby-like vorticity modes and their mutual interaction for arbitrarily 
given radial and latitudinal gradients of the rotation rate. Applying to this equation the 'traditional approximation' 
of geophysics results in a separation into radial- and angular-dependent parts of the physical variables, each of 
which is described by an ODE. The condition for the applicability of the traditional approximation is discussed. 
The angular parts of the eigenfunctions are described by Laplace's tidal equation generalized here to take into 
account differential rotation. From a qualitative analysis of Laplace's tidal equation the sufficient condition for 
the formation of the dynamic shear latitudinal Kelvin-Helmholtz instability (LKHI) is obtained. A small rotation 
gradient causes LKHI of prograde waves (seen in the rotating frame), while strong gradients are responsible for 
retrograde LKHI. The value of the latitudinal rotation gradient has a lower limit, below which LKHI disappears. 
The LKHI result is applied to real solar helioseismology rotation data. It is shown that the m — 1 mode (m = 
azimuthal wave number) instability can develop. This global instability takes place in the whole envelope of the 
Sun, including the greatest part of the tachocline, in radial direction and at almost all latitudes in horizontal 
direction. The exact solutions of Laplace's equation for low frequencies and rigid rotation are obtained. There 
exists only a retrograde wave spectrum in this ideal case. The modes are subdivided into two branches: fast and 
slow modes. The long fast waves carry energy opposite to the rotation direction, while the shorter slow-mode group 
velocity is in the azimuthal plane along the direction of rotation. The eigenfuncions are expressed by Jacobi's 
polynomials which are polynomials of higher order than the Legendre's for spherical harmonics. The solar 22-year 
mode spectrum is calculated. It is shown that the slow 22-year modes are concentrated around the equator, while 
the fast modes are around the poles. The band of latitude where the mode energy is concentrated is narrow, and 
the spatial place of these band depends on the wave numbers (I, m). 
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1. Introduction 

In a recent paper Dzhalilov et al. (2001; paper 1) investi- 
gated which lowest-frequency eigenoscillations can occur 
in the real Sun, moreover, which role they play in redis- 
tributing angular momentum and causing solar activity. 
We found that such waves could only be differential ro- 
tation Rossby-like vorticity modes. However, the general 
nonradial pulsation theory adopted from stellar rotation 
has some difficulities. For slow rotation, when the spheric- 
ity of the star is violated not seriously, the degeneracy 
of the high-frequency spherical p- and g-modes with re- 
spect to the azimuthal number m is abandoned by ro- 
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tation (Unno et al. 1989). Independent of the spherical 
modes non-rotating toroidal flows (called 'trivial' modes 
with a zero frequency) become quasi-toroidal with rota- 
tion (called r-modes with a nonzero frequency; Ledoux 
1951; Papaloizou & Pringle 1978; Provost et al. 1981; 
Smeyers et al. 1981; Wolff 1998). Although rotation aban- 
dons the degeneracy of the modes, it also couples the 
modes with the same azimuthal order, and this makes 
the problem more difficult. For the high-frequency modes 
(e R = lu/2Q, > 1, where u> and f2 are the angular frequen- 
cies of oscillations and of stellar rotation, respectively) this 
difficulty is resolved more or less successfully. For this case 
the small perturbation rotation theory is applied, in which 
the eigenfunctions are represented by power series, the an- 
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gular parts of which are expressed by spherical harmonic 
functions Y" 1 (Unno et al. 1989). These power series are 
well truncated, unless e R < 1, when the role of Coriolis 
force is increasing. 

Namely the low-frequency instabilities are discovered 
in most pulsating stars (Cox 1980; Unno et al. 1989). 
Rotation couples strongly together the high-order g, the 
convective g, and the r- modes with e R < 1 and with the 
same to, but different I (Lee & Saio 1986). Generally the 
matrix of the coupling coefficients to be determined is sin- 
gular (e.g. Townsend 1997). In all papers on the eigen- 
value problem of nonradially pulsating stars, there exists 
a 'truncation problem' for the serial eigenfunctions, the 
angular parts of which are represented by spherical har- 
monics (e.g. Lee & Saio 1997; Clement 1998). 

The governing partial differential equations (PDEs) of 
the eigenoscillations of rotating stars are complicated from 
the point of view of the mathematical treatment, even if 
the motions are adiabatic. This difficulty arises because in 
spherical geometry an eigenvalue problem with a singular 
boundary condition has to be solved. These equations are 
simplified considerably to neglect the tangential compo- 
nents of the angular velocity $7 in the low-frequency case 
£ R < 1 (this means that the motion caused by the Coriolis 
force is primarily horizontal) . This limitation widely used 
in geophysical hydrodynamics (e.g. Eckart 1960) is called 
'traditional approximation' and has been used first by 
Laplace (1778) to study tidal waves (Lindzen & Chapman 
1969). Laplace's equation (or the traditional approxima- 
tion) for e R < 1 is applicable to the stellar case too. The 
main advantage of this approximation is that it decom- 
poses the initial system of equations into a pair of ordinary 
differential equations (ODEs) (e.g. Lindzen & Chapman 
1969; Berthomieu et al. 1978; Bildsten et al. 1996; Lee & 
Saio 1997). The angular parts of the eigenfunctions are de- 
scribed by Laplace's tidal equation. Solving this equation 
numerically by using a relaxation method, Lee and Saio 
(1997) first avoided the representation of the solutions by 
Y^^osO) functions for the =const case, and they had 
no problem with the truncation of the series. 

In the present work for the non-magnetic and non- 
convective cases we receive one PDE in spherical geometry 
for the adiabatic pressure oscillations in the differentially 
rotating star (J7 = f2(r, 8)) with arbitrary spatial gradients 
of rotation (Sect. 2). This general equation is split into 
the 8- and r-component ODEs, if the traditional approxi- 
mation is applied (Sect. 3). The ^-component equation is 
Laplace's tidal equation generalized for the differentially 
rotating case. In Sect. 4 we analyse more qualitatively this 
equation. We find the general condition for the shear in- 
stability due to differential rotation in latitude. We find 
that the smallest rotation gradient is responsible for the 
prograde (seen in the rotating frame) vorticity wave insta- 
bility, while a stronger gradient causes the retrograde wave 
instability. For solar data (small rotation gradients) the 
to = 1 prograde mode instability is possible (Sect. 4.4). 
The possible existence of such a global horizontal shear 
instability on the Sun has been investigated by Watson 



(1981) and Gilman & Fox (1997), that of shear and other 
dynamic instabilities and of thermal-type instabilities in 
stars as well by Knobloch & Spruit (1982) and others. 
Laplace's tidal equation for low frequencies in the rigid- 
rotation case is investigated in detail in Sect. 5. It is shown 
that the eigenfunctions are defined by Jacobi's polynomi- 
als which are of higher order than the Legendre's. 

2. Basic equations 

The fluid motion in a self-gravitating star, neglecting a 
magnetic field and viscosity, may be described in an in- 
ertia! frame by the hydrodynamic equations. These equa- 
tions in conventional definition are written as 



^ + V-( P V)^0, (1) 

9 {jt + V ' V ) V = ~ VP ~ PV$ ' (2) 

P T (i + V -v) a = P^n-V-F, (3) 

V 2 $ = AttG p. (4) 

2.1. Equilibrium state 



We suppose that the equilibrium state (variables with zero 
indices) of the star is stationary and that its differential 
rotation is axially symmetric: 

V (r,8)=flxr = Clrsinde^, (5) 

where fl(r, 6) = fi r e r + He eg, with the components 
fi r = £!cos(9, fig = — rising, and £1$ = 0, is the stellar 
angular velocity of rotation described in spherical polar 
coordinates, (r,8,<fi). Here e.; with i = r, 8, or 4> are the 
unit vectors. We will not include convective motion and 
meridional flows into the initial steady state. In that case 
we may obtain, in particular from the Eq.(|^) of motion, 
the hydrostatic equilibrium relation 

Yll = y [ $ - I|n X r\A + Qr 2 sm 2 8Vn = 0.(6) 
Po V 2 / 

It follows that the effective gravity g cannot be a potential 
field if differential rotation Vf2 ^ is present. This is im- 
portant for rapidly rotating stars where the configuration 
is deformed by the centrifugal force as well as by differ- 
ential rotation. For slowly rotating stars (the Sun as well) 
we may assume that the initial state is only marginally 
disturbed by rotation and g w g = V<I>o can be applied. 
That is, non-sphericity is not essential for the generation 
of waves (Unno et al. 1989). 

2.2. Equations of oscillation 

Small amplitude deviations from the basic state of the 
star may be investigated by linearizing Eqs. (Q-|l). For 
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Eulerian perturbations (variables with a prime) the equa- 
tion of motion becomes (Unno et al. 1989) 







ei—Vi + 2Ux V + e<f,r sin0(V • Vfi) 



Po Po PO 



Here the operator 
Of ~0~t + 06 



(7) 



(8) 



represents the temporal derivative referring to a local 
frame rotating with an angular velocity £1 = £l(r,9). For 
low- frequency waves Saio (1982) has shown numerically in 
detail, that the Cowling approximation is good enough in 
most cases. Thus we will neglect perturbation of the grav- 
itational potential in Eq. (0), $' = 0. We are interested 
in very slow motions such that v p h *C c s , where v p h is 
the phase velocity of the waves and c s is the sound speed. 
Then the incompressible fluid motion limit, c 2 — > oo (it is 
within the adiabatic approximation) may be applied, and 
instead of the Eq. (|l|) of mass conservation we use 



V- V 



1 d{r 2 V r ) 1 0(V e sin9) 



Or 



r sin ( 
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1 0V A 



0<p 



r sin ( 



= 0.(9) 



We have shown in Paper 1, that nonadiabatic effects are 
of great importance for the dynamics of low-frequency ro- 
tation modes. However, here we shall restrict ourselves to 
the adiabatic case only because the mathematical treat- 
ment of wave equations in spherical geometry is rather 
difficult. For adiabatic waves we receive from Eq. (||) 



Op' 



N 2 



1 Op' 



g^-pO — V r -^—, 



Of 



where the squared Brunt- Vaisala frequency 

1 1 dp _j_ dpo 

po dr 



N 



1 \T 1 po dr 



(10) 



(11) 



is written for the slow rotation case where g *=s g. For rapid 
rotation N 2 = N 2 (r, 0) and g should be changed here to 
g. In the incompressible limit (c 2 — > oo) Eq. (Mm reads 



0£ 
Of 



Po — Vy. 



(12) 



Thus we have a complete set of Eqs. ([?], |[ [TJ]) to describe 
adiabatic, low-frequency, non-radial oscillations in a dif- 
ferentially rotating star. For our axisymmetric stationary 
initial state we may represent all the perturbed variables 
V, p' , and p' in the inertial frame as 



V(r, 6, <f>; t) ==> V (r, 0) e i ( m *- h ">*> _ 



(13) 



Considering 0/0t' = — iuj and 0/0<j> — im we find from 
Eq. (|sj) the relation between the frequencies in the inertial 
and rotating frames 



u = luq — fn ri (r, i 



(14) 



from where we get Vw = — mVfi. If we separate the vari- 
able part of the rotation frequency, 0(r, 0) = CIq + Cl(r, 0), 
then oj = luo — mf2o — mQ = ojq — mf2 = luq — fc^fph 
{k$ and w p h are the local azimuthal wave number and 
the phase velocity, respectively). We will study the case 
Cl < fi . From now Q w Q and Vfi = Vfi ^ will 
be used. Low frequencies seen in the rotating frame mean 
that we are close to the resonant frequencies in the inertial 
frame (w « k^Vpb). 

Now excluding p' from Eq. (0) by using Eq. (|l2|), taking 
the projection of this equation onto the rotation axis f2 
and two tangential components as well, and adding Eq. ^ 
we get 



1 5- I V r cos — j Vg sin ( 



1 / d sin 9 

— — [ COS0 ^--J ™ )P 

poiui \ or r 00/ 

2ficos0 T , 1 Op' 

Vg + ; V* = 

IUJ rpoiu 00 

( 0Q\ 
luVj, + j sin I 2tt + r— 1 V r + 



„^ „ . „0fl\ im 

2ficos6» + sin6l— Vg = 

09 J rpo sin t 



V 



imV + S M |_ (r v r) + ±( sin ev e ) = 0. 

r Or 09 



(15) 
(16) 

(17) 
(18) 



Here we have introduced the special parameter j to switch 
to the traditional approximation (Tig — 0, fl r ^ 0). In the 
general case j = 1, and for switching to the traditional 
approximation we put j = 0. Further we will obtain one 
additional equation for p' . 

Let p = sin 9 be a new independent variable and 



£r = 


2ft 


, & 


a = 




-(1 


Oi = 


1 + 


J 



2n Or 1 ' »~ 
(1 - M 2 )(1 + /3 AI ), 

1 + 



a± = p 



£r 

2 1 + A 



rn'op' 



a 

a 3 = -2, 

fc R 



(19) 



,05 = — + 

a e R 



e R is the Rossby number; we are interested in e R < 1. 
Using 



iVj, = jaipV r + a 2 Vg cos ( 



m p 



PoVjJ rp 



(20) 



we get three equations from the Eqs. (|15|-|18|) for V r , Vg, 
and p' . From those Vg is excluded by 



Ve . T> 1 ( 

= jaiV r \a 5 - p— I P, 

cos 9 0,3 \ Op, 

where 

p' = prpoiuP , V r = pV r , 



(21) 



(22) 
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and we get two equations for V r and P: 



N 2 \ ~ 

1 ~ - ja 4 VV = 

d ■ 1 \ d 

0\ + r— - ja 6 - 3 1 li- 
ar \ a 3 J dp 



P. 



d . ( d 
+ r ~dr +3 [ a8 + a9yU a^ 



Vr = A(fl)P, 



(23) 
(24) 



A( M ) = 



1 



a 7 + (l-^ 2 )/i 



_5_ 

dp 



1 / a 

— a 5 -p 

03 V 



(9// 



Here the dimensionless coefficients arc defined as 



0. G 



m 



me R + (l-M i )(l+/3 /i )( — ^-1) 



£r 

J', , 1 n..2 



a 7 = — (1 + /3 M )(1- M 2 ) + 1-2 M 
£r 



08 



; a g = -(l- M 2 )(l + /3 r ), 



(25) 



a* = (1 + /3 r ) [a(3 - V 2 + me R ) - 

-(1- M 2 )M - 2a(l - /z 2 )^ 
, , . m r dp 

bi = 1 + S2 p /3 r , Sip = — . 

£r Po dr 



(small j3 r and (3S) these examinations are also valid. Thus 
for j — Eq. (pq) becomes: 



a 2 /i 2 



a£ R (l - /Li 2 )//" 



P = 



= -I (r 2 ^ 
ip \ dr 2 

with the parameters 

2-3p 2 



d_ 

dr 



1" — - + bir— +b 5 ) p' 



2 - p 2 - e\ 



1 — /i 2 



+ /3 M [2/i 2 - I - 2(1 - /z 2 )(l + (m/e R + 1) 

02 = fro + fti^ + h 2 /3l, 

h = e R (m 2 - l)(e 2 K + p 2 - 1) - 

- M 2 [m(l-M 2 )+4(m-2e R )], 

hi = (1 - M 2 ) [-2m/i 2 + e R (m 2 - 4 + 5/i 2 ) + 

+ 4me R ] — m£j, 

fa = (l-M 2 ) 2 (3m-2e R + 2m/3 M )- 

- (l-//>(l-2rae R ), 



&2 
63 



Deriving these equations the second derivatives CI" (6) and 
n" (r) have been omitted as very small quantities. Eq. ( p3| ) 
allows to obtain one equation for p' : 

ipi$ 2 r + ^2&l + i>A + + faKh + W| p = 0. (26) V> 



A 



„ m 1 m 
= 2 + 3ae p -/3 r — 3 + /3 r — 
£r V £ R 

iV 2 P , N 2 m n 

= 1 - — , S r ^b = -2— — /3 r 



w 2 £r 



r oW 2 
w 2 dr 



(28) 



(29) 



This is the main singular PDE for nonradial rotation- 
gravity waves in a differentially rotating star. The coeffi- 
cients ^1-6 of Eq. (|26| ) are rather complicated, we present 
them in Appendix A. The operators are defined as 







d 



$r = r-sr, 5u = A*^— , and <5 2 = S r S r . 
or dp 



(27) 



First we will study this equation in different simplified 
approximations. The most popular case is the traditional 
appraximation. 

3. The 'traditional approximation': j = 

Strictly speaking, the condition j = in Eq. ( p6| ) is not 
applicable in two points: (1) at w 2 = N 2 (r), that is in 
the turning points in radial direction; (2) at a — or 
e R w cos 2 0, that is in the latitudinal turning point. The 
last one is more important, because the traditional ap- 
proximation niters out such important and interesting 
phenomena as the trapping of Rossby-like waves around 
the equator. In geophysics this phenomenon is investigated 
separately (Pedlosky 1982; Gill 1982). The applicability of 
the traditional approximation for rigid rotation has been 
checked by numerical modelling as well as by experimen- 
tal verification. For minor differential rotation of the star 



93 = £r(1 - V 2 )[a{l - 2a ) - qi], a = 1 - /3 M — , 



£r 



04 = £r(1 - l-i 2 )[a{al - S^a ) + a qi] - q 2 , 



„ m / „ m 
-0 M — 1 + j9„— 



I-261+62, 6 r b 1 = ee p -l3 r — [1 + p r — 

£r V £r 



65 = bf — S r bi - 6i6 2 + b 3 . 

Remembering that Cl(r,8) -C 0©, the left-hand side of 
Eq. ( p8|) is a function of /i = sin0, while the right-hand 
side is a funcion of r only. In that way we may separate 
the variables 



p'(r^) = Q( fJ ,)Q(r) 

Now putting Eq. (| 
V-equations: 



into (C 



(30) 

we receive the '0'- and 



a£ R (l-^ 2 )^ 2 -^-j+<73/i— +q 4 + A 



d 

dp, 



2 2 
a pL 



£r 



— + 6 4 r— + 65 - A^ Q(r) = 



dr 



dr 



e(/x)=o, (31) 



(32) 



These two equations are connected to each other by two 
common spectral parameters: w - the oscillation frequency 
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and A - the separation parameter. Both must be searched 
for as a solution of the boundary value problem. So far 
the logarithmic gradients of the rotation rate are arbitrary 
functions, /3 M = /^(/x), r = /3 r (r). 

4. Generalized Laplace's tidal equation 

Eq. ( |3l| ) is the generalized Laplace equation if differential 
rotation is present, (3^ ^ 0. For rigid rotation, /3 M = 0, 
Eq. ( |3l| ) becomes the standard Laplace equation: 



d 2 



i -si 



F 2 



p 2 dp 2 ' {e\ - ill) 2 d/x* 



1 



2 



A 



6 = 0, (33) 



where /i* = cos#. A peculiarity of this equation is the 
presence of three singular points: at the pole, at the equa- 
tor, and between both if e\ < 1 (our case). Therefore it 
is hard to solve such an equation analytically or numer- 
ically to find the eigenvalues. Since the time of Laplace 
in geophysics investigations were focused on almost two- 
dimensional (horizontal) motions in strongly stratified flu- 
ids with V r ~ 0. In this situation the r-component Eq. ( |32| ) 
does not appear, and one is looking for the eigenvalues A 
in Laplace's Eq. ( |33| ) for a given e R (it is expressed through 
the thickness of the fluid layer, e.g., in a shallow water- 
wave system). For the special cases and for the general 
case too, when the eigenfunctions are expressed through 
the Hough functions (essentially these are the same infi- 
nite series of Y™ harmonics) references could be found in 
a paper by Lindzen & Chapman (1969). 

In astrophysics the r-component Eq. ( [32] ) appears, and 
two equations must be solved together to find both spec- 
tral parameters. In the rigid-rotation case Lee and Saio 
(1997) looked for A numerically in an approach similar to 
that in geophysics, fixing e R in Eq. |3^). We offer here an- 
other approach, where we will find A from the r-equation 
for a given e R < 1. 

It is convenient to introduce into Eq. ( |3l| ) the new 
variable x = /j 2 = sin 2 9: 



where 



d 2 Q 
dx 2 



de 

dx 



A 2 e=o 1 



(34) 



A x = (l-a;)[2-a;-3e 2 i +/3 M (4a;-3-2TOe R )]- 

- 2/3 2 (l-x) 2 + 4, (35) 

77? 777 

A 2 = (l- x )[ m 2 (l-^) + — x + ^ — (4x~3-2el)}- 

- me R (l-a;)+m£ R (l-me R ) + -^-a;(x-a) 2 (l + / 9 AI ) 2 , 



a = 1 - 



This equation determines the tangential structure of the 
eigenfunctions, while Eq. ( |32"| ) is responsible for the ra- 
dial behavior. A detailed investigation of Eq. (|32|) is not 



included in this work. A similar equation for a realisti- 
cally stratified model of the Sun has been investigated 
by Oraevsky & Dzhalilov (1997). We remind that the ra- 
dial structure of the eigenfunctions depends on the sign 
of A (either radiative or convective modes). Now we can 
already estimate an approximate value of A: 



A 

F 2 

R. 



Nlel/u 2 NfjAtl 2 



n 2 x 1CT 6 , 



(37) 



where n is the radial harmonic number, N m is mean value 
of the Brunt- Vaisala, frequency in the radiative interior or 
in the convection zone, and N rn /2tt is the Prandtl number. 
An estimate of Eq. |37]) is done for a solar model, but we 
think similar values of the Prandtl number are valid for 



most other stars too. Thus we can omit from A 2 in Eq. (35 ) 
the last term, if the radial number n is not too large. 

4.1. Fluid velocities 

From Eqs. (|2(], [2l], ^3|) we can derive in the traditional 
approximation the following formulae for the components 
of the fluid velocity: 



V r = 



1 dp' 

Poitoip dr ' 



(38) 



Vg = 



±1 



l-x 



1 



™_2x— )p', (39) 

rpoiuj l+P^a-x y/ x (l-x) V £ R dx - 



1 



rpouj ^Jx 



l-x / d 

2e R x— 

a — x \ ox 



P 



(40) 



Here the different signs ± of Vg correspond to the northern 
and southern hemispheres, so that cos 9 = ±\/l — x. Our 
further aim is to find such solutions for p' — Q(x)Q(r) 
that all the components of the velocity remain limited at 
the pole (x — 0), at the equator (x = 1), and in both 
turning points, where x = a and tp( r ) = 0. 

4.2. Heun's equation 

Now we will impose a restriction to 2/3 M = 
<9(lnf2)/<9(ln/i) w const, the logarithmic latitudinal 
gradient of the rotation frequency. We might take the 
linear dependence /3 ~ x, but for such a profile the 
structure of the solutions is not changed qualitatively. 
Let us introduce the new dependent variable 



6 = x"Y{x), 
where 



2d 

Si = 

So = 



(3^ ± ^(3 2 S 2 - S 2 

5 + 2me R + 2/3 M 
2(1 + ^-4) ' 
^f(3 + 24)-m 2 ( 



2cti: 



(41) 



(42) 



1 + /3„ - el 
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Then for Y(x) we get a new equation from Eq. (p 



x(l - x)(x - a)Y" + 



1 



4<t(1 -x)(x- a)- 



[(a? - lH + e a wi]y = 0, 



(43) 



with 

vq = 
v\ = 



l + 4/? M 
4(1+^) 
m(me R — 1 



H- 2^+0-^(51-1)-^ 
+ 2cte r 



(44) 



4(l + /3„) 

Eq.([l3|) is Hcun's equation (Hcun 1889) in standard form 

x(x - l)(i - a)Y" + [y(x -l)(x-a) + Sx(x - a) + 

+sx(x - 1)] Y' + aft(x - h)Y = 0. (45) 

The Riemannian scheme for this equation is 



P ■ 



1 a oo 
a x 
1-7 1-5 1-e ft 



where the exponents are connected by Riemann's relation 
■ft-j-S-e + l = 0. 



a - 



In the Riemann scheme the first row defines the singular 
points of Heun's equation, while the corresponding ex- 
ponents are placed in the second and third rows. These 
exponents are 



1 - 7 = /? M Si - 2a, 1-5 = 



2 + ft, 1 - 5i 



1 

2' 
3/2 



1 + /3p 



2a 



5+g, 2f3 = S-q, q 



h = l-e a — , 5" = 



■■2a -2- ft x 



a/5 2 -4^0, 
3/2 



(46) 



(47) 



Note that the second exponent at x = a is (1 — e) — > 2 
if e R -> or if ft„ -> 0. If -> oo then (1 - e) -> 7/2. 
That means, the second independent solution of Eq. fl45| ) 
with the exponent (1 — e) is regular at the singular point 
x = a for all variables, which follows from Eqs. (|38|-^o|). 
The exponent (1 — S) = 1/2 also provides limited V$ and 
at the equator (x = 1). The singularity x = oo in the 
Riemann scheme does not occur in our task. The situa- 
tion is more complicated around the pole x = with the 
exponent (1 — 7). Let us consider this point in detail. 

4.3. Condition for the latitudinal Kelvin-Helmholtz 
instability 

If we put the solutions with the exponents and (1 — 7) 



into Eq. (|lj) we get p' ~ 6 ~ x ai - 2 . Then V e ^ ~ x 



(2<7-l)/2 

means that for the regularity of the solutions at the pole 
the condition Re(2a) > 1 must be obeyed. On the other 
hand, an instability is possible when the eigenfrequencies 



are complex, that means complex a. For the latter it fol- 
lows from Eq. that the necessary condition is S2 > 0. 
It is clear that the axially-symmetric mode with m = is 
excluded. For lower values of the rotation gradient |/3 M | < 1 
the necessary condition S2 > demands for the prograde 
waves (me R > 0) the condition ft, > 0, which is more 
realistic for stellar situations (equatorward spinning up at 
the surface with radius r). Rayleigh's necessary condition 
for instability (Rayleigh 1880; Watson 1981) says that the 
function Rl = d 2 (11 sin 2 8)/dcos 2 8 (gradient of vorticity) 
must change its sign in the flow. Rewriting this function 
in our definitions we get that 



Rl 



M 2 (l + 4/^)] 



(48) 



may change its sign if ft^ > 0. There are instability pos- 
sibilites for negative ft, which are not considered in this 
work. However, all formulas are valid for this case too. 
The sufficient condition for instability is obtained from 



Eq. (42) and reads ft 2 S 2 < S%. The regularity condition 
at the pole ft^Sx > 1 can be rewritten as 



< 



f 2 

fc R 



ft^ftp + 3/2) 



X3- 



(49) 



By this condition the phase space {e R /m, e R } is divided 
into three parts, depending on the values of ft For 

< ftp < 1/2 we have the following situation: if £ R < 

1 - ftfi(ftn + 3/2) the condition Eq. @ is fulfilled for pro- 
grade waves e R /m > (region I); in the opposite case 
when e 2 > 1 - ft^ft^ + 3/2) the condition @ is fulfilled 
for retrograde waves with Er/to < —1/2 (region III); for 
e 2 , = 1 — ftfi(ftfi + 3/2) these regions are separated by the 
asymptote \3 = 00 ■ 



For ftp > 1/2 (strong gradients) the condition Eq. (|49| ) 
is met only for retrograde waves in the range —1/2 < 
£r/to < (region II). For ft, = 1/2 we get for any e R 
that \3 — — 1/2- This is the line between the regions II 
and III. All three regions are shown in Fig. la. It is seen 
that the regularity condition is working for |m| > 1 and 
|£r| < 1 if ftn 7^ 0. For very small ft^ only modes with large 
m are possible. The smallest m ~ 1 modes may appear in 
the limit e R wl. These conclusions are correct only if the 
instability occurs. 

Now let us consider the second condition, the complex 



frequency condition ft,Sf 
rewritten as 



< 52- This inequality may be 



X2 < 



Xl,2 



£r 

— < Xi, 
m 



(50) 



^(5 + 2ft,) 2 



&* + 



b 2 - a*e R (5 



2/3, 



(1 
(3- 



4)(i- 
24) a 



■0u 



[ft, (5 + 2ft, 



In the limiting cases we have 



Xi 



12(1 + ^) 



X2 ~ for e R 



(51) 
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1.0 



0.8 



m 




■ o.wz 



0.1 0..05 



Prograde waves. 




Fig. 1. The domains of validity of the solution regular- 
ity and of the instability conditions in the phase space 
{e R /m, e R } for given values of the rotation gradient (3^. 
(a) shows the values e R /m = % 3 , Eq. (^|). In the area 
Er/ti < X3 the solutions are limited at the pole, (b) shows 
the behavior of xi (dashes) and X2 (solid). Between the 
solid and dashed lines with the same labels (values of /3 M ) 
LKHI is possible. 



Xi 



0, X2 



8/3^ 



(5 + 2/3 M ) 2 



for £ R 1. 



(52) 



In Fig. lb we plot for \m\ > 1 and e R < 1 the curves \i 
and X2 versus e R for a wide range of /3 M . A comparison of 
Figs, la and lb shows that in the region III with e R /m < 
— 1/2 LKHI will never appear. In the region I LKHI is 
possible for prograde waves, if /3 M < 1/2 and e R < 1 — 
PuiPn + 3/2). With decreasing [3^ the solid and dashed 
curves for the same (3^ are close to e n ~ 1- For retrograde 
waves LKHI is possible for strong gradients of [3^ > 1/2 
only in the range —0.2 < £ R /m < 0. Here (x2)min = —0.2 




-0.20 
0.0 



Fig. 2. The wave instability areas (hatched) in the phase 
space from an overlap of Figs, la and lb. The labels in 
the areas are the /3 M values. Both prograde and retrograde 
waves (3^ are limited: for prograde waves the instability is 
possible if (3^ < 0.5, for retrograde waves the instability is 



possible if (3^ > 3.5. 



is valid for e| 



1 and [3^ — 5/2, which follows from 
Eq. (|2j). 

The total condition for the existence of spatially stable 
but temporarily unstable waves reads as follows: 



X2< — < min(xi,X3)- 
m 



(53) 



Figs. 2a, b show the validity ranges of this condition 
for some typical values of (3^ . The hatched areas are places 
where LKHI is possible. These figures are obtained by 
overlapping the Figs, la and lb. For prograde waves on 
both sides these hatched areas become very narrow: with 
decreasing (3^ the extent of the hatched area decreases and 
tends to the point {e R , m} = {1,1}. We will see that this is 
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the solar case. The hatched areas of LKHI disappear with 
decreasing /3 M . This means we have a lower limit (/3 M ) m i n - 

For retrograde waves it is sufficient to write the condi- 
tion as xi < £ R./ m < X3- I n Fig- 2a xi is the dashed curve 
and X3 is the solid curve. LKHI is possible if /? M > 3.5 for 
> £ R /rn > -0.2. 

In Fig. 2 the hatched areas means that only these 
hatched areas for given /3 M are possible, if LKHI takes 
place. Outside these hatched areas regular solutions are 
impossible. The case without LKHI (neutral oscillations) 
must be investigated separately. 

4.4. Solar rotation profile 

Let us consider at which places we might expect LKHI in 
the Sun. Unfortunately, it is not clear how the core rotates. 
Nevertheless some rotation gradients close to Sun's center 
might exist, and we could expect LKHI there. It is known 
from helioseismology that the radiative interior has a very 
small (3^, but the exact value is unknown. We have better 
information on the rotation profile of the solar envelope, 
including the tachocline. Helioseismology data may be de- 
scribed by different approximate formulae. One of these is 
(Charbonneau et al. 1998) 



Q(r,0) 

n, 

tt c /2lT 
Cl 



n c + i[i + erf(A)] (n,(0)-n c ), 



(54) 



a. 



ci cos 9 + c 2 cos 8, A = 



432.8 nHz, ft oq /27r = 460.7 nHz, 
-62.69 nHz, c 2 = -67.13 nHz, 



where r c = O.713i?0 is the radius at the bottom of the con- 
vective zone, and w = 0.025Rq is the tachocline thickness. 
We can easily check that our approximation £l/Q c <C 1 
is always applicable, if Eq. ( |54| ) is represented by Q = 
fl c + f2(r, 6). The maximum value (f2/f2 c ) max ~ 0.06 for 
= it/2 (equator) is in the convection zone. From Eq. ([54]) 
we find the latitudinal gradient of rotation 



0u = 



1 + erf (A) 
2tt 



sin 9(ci + 2c 2 cos 



(55) 



Our supposition about /3 M ~ const and (3 r ~ const is 
based on the presentation ft £li(r) + SI2 ■ We could 
receive such an approximate formula instead of Eq. (pi}) , 
but we need only local values of the gradients for which 
Eq. d55| ) is acceptable. Using this formula we show in Fig. 3 
the (3^{6) dependence for different r/R@. We see that in 
the Sun [3^ < 0.02, and the maximum is in the photo- 
sphere. From Fig. 2 follows that the LKHI of retrograde 
waves is not present in the solar case. LKHI of prograde 
waves in the Sun occurs in the upper right corner of Fig. 2a 
which is enlarged in Fig. 4. Here we see that the LKHI 
area disappears when (3^ < 3 10~ 4 . This boundary is lo- 
cated at the bold horizont in Fig. 3. Thus the prograde 
waves become unstable in the Sun at those places where 
3 10" 4 < /3 M < 2 10~ 2 . It means that LKHI is possible 
in the area r/i?© > 0.6725 which includes the greatest 




T 



0.69 



0.675 



40 60 80 

Latitude, deg. 

Fig. 3. The local estimate of the logarithmic gradient of 
the solar rotation frequency (3^ for real solar data from 
helioseismology depending on the co-latitude and the ra- 
dial distance (the labels are values of r/R Q ). The bold 
horizont is /3 M = 3 10 -4 , above which the prograde waves 
become unstable (see next picture). 



m 



.3 



2 .1 



0.996 0.998 
2 



Fig. 4. Enlarged part of Fig. 2a for the smallest gra- 
dients of rotation. The labels 1, 2,..., 6 correspond to 
/3 M = 310- 4 ,510- 4 , 110- 3 ,210- 3 ,310- 3 ,410- 3 . Areas 
of instability exist only if /3„ > 3 10 -4 . 



part of the tachocline, the convective zone, and the pho- 
tosphere. With increasing r the LKHI zone expands from 
middle to high latitudes. Figs. 2a and 4 show that LKHI 
is occurs at high frequencies (e R ~ 1) and in global scales 
(e R /m > 0.5). Considering that m 7^ is an integer we 
get to = 1. 

However, our quantitative analysis of LKHI is based 
on the general Riemann scheme of Hein's equation, which 
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is valid only if the middle singular point x = a is far from 
the other edges at x — and x = 1. Thus, the limiting 
cases e R — * and e H — * 1 (the latter is more important 
for solar LKHI) should be considered separately. In these 
limiting cases the regularity condition Eq. (48) may be 
changed, and the curve in Fig. 4 limiting the instability 
areas from below may be shifted. In this case LKHI with 
higher m-modes should be expected. 

5. The low-frequency waves 

After the qualitative analysis of Heun's Eq. (^) we can 
start a quantitive analysis. Note that the qualitative 
conclusions drawn above are valid for the more general 
Eq. ( [34] ) with A term. Heun's equation with four singular- 
ities in the general case is solved by a series of hyperge- 
ometric Gauss functions. A similar task has been consid- 
ered for the damping of MHD waves at resonance levels 
by Dzhalilov & Zhugzhda (1990). We will start to study 
Eq. (Eq) for some simple limiting cases. At high frequen- 
cies (ej ~ 1, when LKHI is acting in the Sun) and at low 
frequencies (e^ <C 1, when the waves are stable against 
LKHI in the Sun) Heun's equation is strongly simplified. 
In these cases the singular level x = a is shifted either to 
the pole or to the equator. For both cases solutions are 
expressed by one hypergeometric function. 

In the present work we consider particularly the second 
case. Let e\ <C 1 + (3^. Then we have a « 1 and h w 1. 
Eq. 



is now the hypergeometric equation: 



x(l - x)Y" + [7 - (a + + l)x]Y' ~ &f3Y = 0, 



(56) 



where all parameters are defined by Eqs. ^?j). In 

these definitions £p should set to zero. Then the 8-part of 
the pressure perturbations, Eq. (41), is expressed by two 
Gaussian hypergeometric functions: 



6 = dx 171 Yi(x) + C 2 x a2 Y 2 (x). 



Here Y 1 .2 (x) = F(a, (3; 7; x)\ a=c 



(57) 



(58) 



and C1.2 are arbitrary constants. This general solution 
includes LKHI for larger /3 M too. This could be realized 
perhaps in other, younger stars. For the Sun we have (3^ <C 
1. We will finish this paper by considering in detail the 
more popular case when rotation is uniform, (3^ = 0. 

5.1. Rigid rotation case 

Using the conditions JL — and e 2 , <C 1 the parameters 
in the solution Eq. ( p7| ) are strongly simplified. Because 
2<7 = ±|m| only a regular solution at the pole (x = 0) 
will be left. In the standard definitions of hypergeometric 
functions (Abramowitz & Stegun 1984) we have 



9 = Cx^/ 2 F(a,b;c;x) 
a= \{-\ 



(59) 



+ \m\ + q I , b 



N 

c 

0.) 




Fig. 5. The spectrum uj(l,m) of low-frequency retrograde 
modes. The numbers at the curves are the / values. For the 
calculation of the spectrum, Eq. (jF" ' 
nHz is used. 



n = n cq = 460.7 



c = 1 + \m\, q 



, and C = const. 

4 £ R 



Note that c — a — 6 = 3/2. The analytical continuation of 
the solution Eq. (|59"|) to the equator, x —* 1, gives 



e = Cx lm \ /2 



AF(a,b;-~;l-x)- 



+ B(l- xf /2 F{c~ a, c -&;-;!- x) 



(60) 



Here the continuation coefficients are (Abramowitz & 
Stegun 1984) 



A 
B 



T(c)r(3/2) 
T{c-a)T{c-b) 

r(c)r(-3/2) 

T{a)T{b) ■ 



F(a,b; c; 1), 



(61) 
(62) 



Eqs. (^9[ p0|) mean that pressure is limited everywhere in 
the hemisphere. Now let us consider the velocity compo- 
nents. Putting Eq. @ into Eq. @ we get 



± 



C x 



r \m\/2 



rpniuj y/x(l - x) 
2ab 



£r 



m\ I F(a, 6; c; x)- 



x F(a + 1, b+ 1; c+ l;ar) 

c 



(63) 



Taking into account that F(a, b; c; 0) = 1, we receive from 
the regularity of V$ at x = that |m| > 1. Axially- 
symmetric waves m = cannot be formed. The contin- 
uation of the solution Eq. (|6|) to the equator, x — ► 1, 
gives 



Vo 



± 



C Q £ R x (\m\-l)/2 
rp iuj 



A 



(64) 
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o.oi 



" 0.00 
_0 

o 

> 

3-0.01 

o 
^_ 




I m I 

Fig. 6. The normalized group velocity v g /2£lr s'm9 as a 
function of the azimuthal numbers for given degrees Z. 



Til \ I 

A = [--|m|jF(o,i;--;l-a;) 

1 



-x4abF(a + 1, 6+ 1; -; 1 - x) 



Co = 



1 ml (l-x)+3x 



F(c-a,c-b;-;l-x) 



4 7 

+ -x(l—x)(c—a)(c—b)F (c—a+ l,c-6+l;-; l-x). 

5 2 

As for x = 1 the functions C\ and C2 are limited, we 
have the relation ^4 = from the regularity requirement. 
Using a property of the gamma-functions (its presentation 
as an infinite product) in Eq. (|^) we get the condition of 
quantization 



21 



m 



I = 0,1,2, 



From here we get the simple dispersion relation 

1 



£r _ 1 IjJ 

m m 2S1 



(2 + 2Z + |m|)(3 + 2Z + |m|) 



(65) 



(66) 



Since £ R /m < for any Z > 0, only retrograde modes (as 
seen in the rotating frame) are possible. 

5.1.1. Spectrum of retrograde modes 

The new dispersion relation Eq. (|6^) completely differs 
from the dispersion relation of the almost toroidal r- 
modes. Their dispersion relation can be derived from 
Eq. @ if we formally set l s = 2(1 + I) + |m|. Then 



-2Q 



i s (i + i s y 



(67) 



However, here the degrees l s are functions of m. Due to the 
coupling of the modes in our case the eigenfuncions can 



never be expressed by the associated Legendre functions 
Pf 1 . For the r- modes the axi-symmetric modes (m = 0) 
with uj = are possible, while our case rejects this case. 
Unlike the r-modes the spectrum Eq. (|6^) has a max- 
imum at m 2 = ml = 2(1 + Z)(3 + 21). If / > 1 we 
have |£R,| m ax ~ 1/8Z. The spectrum of retrograde waves 
is shown in Fig. 5. Rossby waves in geophysics (Pedlosky 
1982) have a similar spectrum. Using Eq. ( p^ ) we may de- 
fine the local phase and group velocities (at fixed latitude 
#0 and radial distance tq) in the azimuthal plane 



20rg sin 9q 



(2 + 2Z + |m|)(3 + 2Z + |m|)' 
m 2 -2(1 + Z)(3 - 



2Qro sin 9q 



21) 



(2 + 2Z+H) 2 (3 + 2Z + |m|) 2 



(68) 
(69) 



We see that the group velocity changes its sign at the max- 
imum of the spectrum, where m 2 = m 2 ,. Fig. 6 shows that 
the I = mode has a maximum group velocity. Long waves 
carry energy opposite to the rotation direction, while a 
packet of short waves is carried in rotation direction. The 
facts that the frequencies oj(m) have a maximum (two dif- 
ferent \m\ correspond to the same uj) and at the maximum 
of uj = cj max the group velocity changes the direction hints 
at the existense of two branches of oscillations. Solving 
Eq. ( |66| ) for \m\ we get these branches. Let ui > and 
m = — \m\. Then 



1 



(70) 



\m\ = - — (wi T \fwi) = mi, 2 , 

wi = 1 - £ R - 4e R (l + 0, 
w 2 = (l-e R ) 2 -8£ R (l+Z). 

From here we get an upper limit for I if e R is given: 
Z < Z max = (1 — e R ) 2 /8e R — 1. For such degrees of Z the 
azimuthal numbers are also limited: mi < \m\ < ni2- 
For £ R — * we have Z max — * 1/8£ R , 1712 —>■ 1/£ R , and 
mi — > 0. However, considering the regularity of the solu- 
tions Eq. ([33j), we must take m\ > 1. From here we get 
the lower limit of Z: 



I > L 



1 



4 

£ R 



7 



(71) 



Since Z > 0, we get £ R < 1/12 from Eq. @. For £ R = 1/12 
we have Z m ; n = 0. Thus the eigenmodes exist for £ R < 1/12 
if Z m i rl < I < Z max and mi < |m| < m 2 . For I = Z max we 
have mi = m2, and mi = 1 for Z = Z m j n . This situation 
is shown in Fig. 7 for different values of £ R . A decrease of 
the frequency decreases the domain of the existence of the 
modes. 

<M) into Eqs. (p 



Setting Eq. 



69) gives 



W\ ± Jvj~2 



251ro sin 9a 

v ± 

257ro sin 9q 



4(1 + I)(3 

W2 ± Wi 



21)' 

W2 



4(1 +Z)(3 + 2Z)' 



(72) 
(73) 



Here v^ h is the phase velocity of the fast modes and 

±-0 



v ph th a t °f the slow modes. For I 



Zmax w e have v„ 
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Fig. 7. The possible domain of the existence of eigen- 
modes for given Rossby numbers e R (the numbers at the 
curves). For example, the case e R = 0.002 (close to the 
22-year modes) is emphasized. All possible values of the 
azimuthal numbers (to = mi and m = mi) are on this 
curve for given discrete I in the range / m i n < I < Z max (see 
text). If e R > 1/12, the eigenmodes disappear. 



,-2 



O 


> 



03 

1, 



,-4 




0.002, - ' 



0.0015 - - 



0.001 



Fig. 8. The phase velocities of fast (solid) and slow 
(dashed curves) modes versus I. The numbers at the curves 
are e R values, similar to Fig. 7. Here the velocities are nor- 
malized to (— 2f2r sin0o)- Every curve is restricted to the 
range / m ; n < I < i max - 



and v+ h 



v ph . In Eq. ( |70P and in Fig. 7 the TOi branch 



corresponds to the fast mode, but mi to the slow modes, 
since TO2 > mi. In Fig. 8 the normalized phase velocities 
(with inverse sign) for the selected values of e R in Fig. 7 
are shown versus /. Both branches are retrograde modes 
(Vp h < 0). Using Oi? ~ 2 km/s for the Sun, we get from 
Fig. 8 very slow phase velocities. The fast wave velocity 



(solid lines) depends more strongly on I. With increasing 
e R both branches are accelerated. 

In Fig. 9 the group velocities are presented in the same 
way. For fast waves the group velocity is always parallel to 
the phase velocity (u+ < 0), while for the slow waves we 
have the opposite behavior v~ v > 0. Slow modes packets 
carry off energy in the rotation direction. Always > 
\v~ T \ is valid. With deceasing e R the range [i m i„,i ma x] is 
shifted to the right-hand side, and it is seen in Fig. 9 that 
v~ for such low s-a is almost zero. 

gr in. 

Note that m = I modes are always fast modes. 
5.1.2. The eigenfunctions 



Taking into account the quantization condition Eq. (65) 
in the solutions Eqs. (60, |64|, and ^) we obtain the 
eigenfunctions. Turning from complex velocities into the 
real displacements, V — —iui£ (recall that V is the veloc- 
ity seen in the rotating frame), we get 



p' = Q(r)Q(9) cos m (<p 
£e = ^cos 



TO I 



U = C sin m ( 

£r = &COS 



U-2—Qt 

V TO 
TO 



U-2—Qt 

V TO 



Here the amplitude functions are 
1 6 dQ{r) 



c = 

Co = 



lu 2 — N 2 po dr 
Q(r) 



4il 2 rp 
Q(r) 



Me). 



(74) 
(75) 
(76) 
(77) 

(78) 
(79) 
(80) 



An 2 rp ' 

The amplitudes are expressed by Jacobi's polynomials: 
B(sin 6>) |m| cos 3 6*^(60, (81) 

(82) 



e 



B(sin< 



)|m|-i cos 2 ( 



TO 4/fc . 2 

— Fi sin 0F 2 

£r 5 



R 21k 

Ua = -sin(20)cos 2 0(sin0)l m l J F 2 (0), 

£r 5 



Fx=F 



-I, k: —: cos 9 

' ' 2 



(83) 
(84) 



= ( _i)< ^ "f(5/2) P (\m\,m { e) 



% (5/2), j\ 



7 



F 2 =F [ 1-1,1 + k;-:cos 



(85) 



, Q - l)T(7/2) p (i+\™\,5/2) ( e) 

[ ' ra + 5/2) 1 ' 
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Fig. 9. Absolute values of the group velocities of fast 
(solid) and slow (dashed curves) modes normalized to 
2f2r sin O versus I for selected e R . For fast waves u+ < 0, 
for slow modes v~ > 0. For I = l max ■ 



we have v+ = w gr = 0. 



i-i 

E 



(l-Q^i + fc), 

(7/2)iil 



(cos 9) 



2.) 



z 5 ,ui,, /« r (/ + i) 



(86) 



and F2 = if 2 = 0. The eigenfunctions of the retro- 
grade low- frequency modes Eqs. (|74|-|77j) are written for 
the hemisphere < 9 < §. In the other hemisphere we 
have to change cos 9 to | cos 9\ and to change the sign of £g . 
We put C = 1 because all eigenfunctions are multiplied 
by the solution of the r-Eq. (|3^) . Changing the sign of m 
the function ^ keeps its sign because £ r /to < 0. It follows 
from Eqs. (^4, 85) that the eigenfunctions are represented 
by a limited number of trigonometric functions. The num- 
ber of terms is defined by I in the range Z m j n < I < 2 m ax- 
This is qualitatively different from the presentation of the 
rotation-mode eigenfunctions by infinite series of spherical 
harmonics. 

Let us consider some particular cases. 

a) Pole (9 = 0): taking into account BF\ — 1 we have 
9 = p' = = ^ = 0. & = if |m| > 1. For \m\ = 1 
we have Cba = -2(3 + 21) (2 + I). 

b) Equator (6 = §): here F x = F 2 = 1 and p' = & = 

& = 0. 

c) I — case: then B = Fi = 1, F2 = 0, 

9 = (sin(9) |m| cos 3 (9, 

S, BA = -(2+ |m|)(3 + |m|)cos 2 0(sin0)l m l- 1 , 

and £0 = (in our e R <g; 1 limit). 

d) I ^ \m\, \m\ is not large, and 9 ^ 0, | as well. 



The latter case is important because the range Z m i n < 
I < ^max is large for small e R . Remember that for e R = 
1/12 we have Z m ; n = 0. Using the asymptotic formula for 
Jacobi's polynomial (or the hypergeometric function) of 
large degree we have 



9 



\m\\ cos A cos 9 



VttI l\ m Wsin9 
\m\\ cos A 
V^I/ |m| (sin 0) 3 / 2 
14 Imllfc sin A 



(87) 
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-fctan^tan A 



(89) 



These extremely simple asymptotic formulae for the eigen- 
functions might be used in most cases. The formula for the 
case of large I with large |m| also could be found, e.g. in 
the book of Bateman & Erdelyi (1953). 

In Fig. 10 for some typical selected (I, m) pairs the am- 
plitude functions, Eqs. (|3l|-|S3"|) , are shown as function of 9. 
The first row is the pressure &(9) function normalized to 
its maximum. The first I — 1 window represents & for dif- 
ferent I to I (increasing \m\ from left to right). As the eigen- 
functions are multiplied by a sin'" 1 ' 9 factor, the ampli- 
tudes are strongly suppressed around the pole. Increasing 
I m I for a given I shifts the maximums toward the equator. I 
is the surface node number of the 0(0) function. Increasing 
I for given \m\ (the second window of pressure with to = 1 
in Fig. 10) contrarily suppresses the amplitudes around 
the equator, and the maximum is shifted toward the pole. 
/ = \m\ is the equilibrium case. In the third window of 
pressure the balance latitude with a maximum amplitude 
is defined by 9 — 20° for all I = \m\ modes. 

The second and third rows of Fig. 10 are the lati- 
tudinal and azimuthal (£^a) eigenfunction ampli- 
tudes, respectively, normalized to the maximum of £9 a, 
see Eqs. ( y% p0| ). The latitudinal amplitude behavior is 
similar that of the pressure. The azimuthal amplitudes are 
smaller than the latitudinal amplitudes, but with a change 
of I a redistribution of the amplitudes will not take place. 

has practically the same amplitude at all latitudes 
and for all I. 

From Fig. 10 follows that we can expect an interesting 
behavior of the eigenfunction amplitudes, when both I and 
\m\ are large. A suppression from two sides may evoke a 
concentration of wave energy in narrow latitudinal bands. 
For example, this is the case for the 22-year solar mode. 

For the 22-year modes we take W22 = 27r(1.441 nHz), 
for which e R w 0.0016. Then we derive from Eq. ( |7l] ) the 
limiting values of the integer /, 11 < / < 76. For all I in 
this range we find from Eq. ([70]) , rounding off, integer az- 
imuthal numbers toi and mi for the fast and slow modes, 
respectively. Putting these integer numbers into Eq. (66) 
we find the deviation S = LO22 — w from the central fre- 
quency due to the integer azimuthal numbers. The results 
are given in Table 1 in Appendix A. It is seen that the 





-1.0 1 1 -1.0 1 l -to t 

20 40 60 80 20 40 60 80 20 40 60 80 

Latitude (in degrees) Latitude (in degrees) Latitude (in degrees) 

Fig. 10. The amplitudes of eigenfunctions versus latitude for given pairs (I, |m|), Eqs. fl8l])-([83|). The pressure row 
shows the Q(9) function normalized to its own maximum. From left to right all curves are related to the wave numbers 
given above. The middle and last rows are similar to the first row, but show now the latitudinal (£9 a) and azimuthal 
(£,4>a) amplitude functions, normalized to the maximum of (£qa)- 



fast modes with low I have larger deviations. This table 
includes all possible (I, \m\) pairs which correspond to the 
22-year period. For some example pairs of (I, \m\) we plot 
in Fig. 11 the latitude dependence of the the quantity 
(£g +£^) 1 / 2 , averaged over the wave period, which charac- 
terizes the energy density of the modes. The hemisphere 
is divided into two equal parts: slow modes are located 
around the equator (solid lines), fast modes are concen- 
trated around the pole. Each (I, \m\) pair is located in a 
narrow latitudinal band. Note that the slow modes (the 
group velocity of which is in the rotation direction) with 
sunspot-like spatial scales are at latitudes of 30 — 40° from 
the equator. 



5.1.3. Flow patterns 



The eigenfunctions Eqs. (|74 - 77 ) allow us to discuss the 
flow character produced by the waves, even if the solu- 
tion of the radial equation Q(r) is unknown. Excluding 
from these equations the time-dependent phase we can 
receive the trajectory equations of the fluid elements. In 
the meridional (r, 0) plane we have 



§1 = .2 Q(r)ff>(r) 6 



6A 



rQ'(r) 6 



tan(a r ). 



(90) 



It follows from here that the motion in the meridional 
plane is linear. The poloidal displacement vector in the 
meridional plane is inclined to the radius by an angle a r . 
For a r sa J we have meridional flows and for a r ~ 
radial flows. However, since £g/£ r ~ l/sin(20), for any r 
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22-year modes 




20 40 60 

Latitude (in degrees) 

Fig. 11. The normalized energy density of the 22-year fast 
(dashed) and slow (solid) modes versus co-latitude. The 
numbers at the curves are (/, \m\) pairs taken from Table 
1. 

the motions around the pole and around the equator will 
be almost meridional (\t;g\ 3> |£ r |). We see from Eq. ( |90| ) 
also that the direction of the flow vector will be changed to 
7r/2 if we pass along the radius from the node of the Q(r)- 
function to Q'(r). Note that every node of Q'{r) is located 
between the neighbouring two nodes of Q(r). For large 
radial numbers (such orders are expected) these nodes are 
located close to each other. Then we have a complicated 
motion in the meridional plane. 

At the surface of the cone over (r, <f>) the trajectory 
of each fluid element is an ellipse around the equilibrium 
point. The cone displacement equation is 



a 2 



f 2 



= l. 



The ratio of the semiaxes £1 and £* 



q Q{r)^{r) Alk F 2 . 

— = — £r — sm9 

£ rQ'(r) 5 F x 



(91) 



(92) 



means that close to the axis of rotation the flow is directed 
along this axis. On the bigger cone through the equator, 
where F± = F 2 = 1, the motion is mainly azimuthal since 
Ik > 1. An exception are in this case the nodes of the 
Q(r)-function at the radius. For large I ^> 1 the ratio 
~ tan(2i0)/ cos 9 means that the character of the 
flows is changing at middle latitudes. 

At the surface of any sphere with a radius r the motion 
of the fluid elements is on trajectories of the ellipse 



a 2 



ii. 

?c4 



1. 



(93) 



Here the ratio of the semiaxes is independent of the radius: 



Uk 
5e R 



Fo 



Fi 



ilk 



Fo sin 2 9 



(94) 



It is seen that near the equator and near the pole the 
motion is mainly in e# direction. For large I 3> 1 this 
ratio is 



Uk 

5e R 



sin 9 tan A 
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k tan 9 tan A 



(95) 



The flow pattern considered here drifts opposite to the 
sense of rotation with a speed u p h, if it is observed in a 
frame rotating with the star. 

6. Conclusions 

In the present work we have derived the general PDE gov- 
erning non-radial, adiabatic, long-period (with respect to 
the rotation period) , linear oscillations of a slowly and dif- 
ferentially rotating star. This general equation includes all 
the high-order g-modes and all possible hybrids of rotation 
modes as well as their mutual interaction. The geophysi- 
cal 'traditional approximation' considerably simplifies this 
general equation, and we get two ODEs for the r- and 
(^-components instead of one with arbitrary gradients of 
rotation f2(r, 9). We have received a more stringent con- 
dition for the applicability of this approximation to the 
pulsation of stars. Only for very low frequencies this re- 
striction is the same as that of the standard case. 

The ^-equation is Laplace's equation generalized to 
the latitudinal differential rotation. Without solving this 
equation qualitatively we found the exact condition for the 
appearance of a global instability. This instability is driven 
by the latitudinal shear, it is not influenced by buoyancy. 
We call that a 'latitudinal Kelvin-Helmholtz instability' 
(LKHI). The appearance of LKHI strongly depends on the 
Rossby number (the ratio of rotation period and period 
of motion), on the azimuthal wave numbers and on the 
latitudinal rotation gradients. Very large gradients pro- 
duce retrograde waves (seen in the rotating frame), while 
a slower rotation gradient is responsible for prograde mode 
LKHI. The rotation gradient has a lower boundary below 
which LKHI is not possible for any Rossby number or az- 
imuthal number m. 

We have applied the LKHI condition to the hclioscis- 
mological data of the Sun. Here a global LKHI is pos- 
sible for the m = 1 mode at practically all latitudes. 
Radially the LKHI is extended from the greatest part of 
the tachocline up to the photosphere. The LKHI for the 
Sun was first obtained by Watson (1981). According to 
his results the instability is possible only at photosphere 
layers. Later Gilman & Fox (1997) have shown that such 
an instability is possible in the tachocline too, if strong 
toroidal magnetic fields are included. Our results show 
that the instability of the m = 1 modes and other modes 
is possible without magnetic fields, in contradiction to 
Gilman & Fox (1997). This difference is probably con- 
nected with the incompleteness of the equations used by 
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Watson (1981) and by Gilman & Fox (1997); their equa- 
tions are two-dimensional only. 

The exact solutions of Laplace's tidal equation for 
lower frequencies are expressed by Jacobi's polynomials. 
Just for lower frequencies the numerical calculations of 
stellar pulsation analyses meet great problems, when look- 
ing for the cigenfunction as infinite scries of Lcgcndrc func- 
tions. The eigenfunctions, defined by higher-order polyno- 
mials of Jacobi, cannot be expressed by convergent series 
of associated Legendre functions. Every Legendre function 
is a particular case of a Jacobi polynomial. 

It has been shown here that the retrograde (slow and 
fast) modes with high surface wave numbers (l,m) are 
energetically concentrated in narrow bands of latitudes. 
This analysis was done for the 22-year modes as an ex- 
ample. Such a concentration of mode energy in a narrow 
spatial area makes such modes vulnerable to different in- 
stability mechanisms such as the e-mechanism considered 
in Paper 1. 



Appendix A: Coefficients of the main equation 

In Section 2.2 our main Eq. (26) for the pressure pertur- 
bations has been derived: 



[ipiS* + ifaSfr + 1p 3 5 r + V> 4 ^ + ip 5 S r S^ + il> & ]P = 0. (A.l) 
This singular PDE has the following coefficients: 
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where 



fi = ipi(2 + jas) - 5 r ipi - jagSpipi, 

Tfl 

f 2 = l + av f3 r -ja 6 , 

£r 

a* = ( / i 2 -l)(l+/3 AJ )/a=l-l/a 3 . 



Taking into account Vw = — m\7fl we can easily obtain 
all the required derivations of the parameters. 
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Table .1. Results for the 22-year period: frequency deviations 5 = W22 — w for permitted quantum numbers (I, | 
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